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PATENT 

SYSTEM AND METHOD FOR MEASURING 
WAVE DIRECTIONAL SPECTRUM AND WAVE HEIGHT 



5 Background of the Invention 

Field of the Invention 

The present invention relates to underwater acoustic measurement systems and, 
more particularly, to acoustic Doppler current profilers used to measure wave directional 
spectra and surface wave height. 

10 

Description of the Related Technology 

The use of Doppler sonar to measure currents in a fluid medium is well 
established. Conventional acoustic Doppler current profilers (ADCPs) typically use an 
array of acoustic transducers arranged in the well known Janus configuration. This 

15 configuration consists of four acoustic beams, paired in orthogonal planes. The ADCP 
measures the component of velocity projected along the beam axis, averaged over a 
range cell whose beam length is roughly half that of the emitted acoustic pulse. Since 
the mean current is assumed to be horizontally uniform over the beams, its components 
can be recovered simply by differencing opposing beams. This procedure is relatively 

20 insensitive to contamination by vertical currents and/or unknown instrument tilts. 

The analysis of waves in a fluid medium is much more complicated, however. 
Although the wave field is statistically stationary and homogeneous, at any instant of 
time the wave velocity varies across the array and as a result it is not possible to 
separate the measured along-beam velocity into horizontal and vertical components on a 

25 sample-by-sample basis. If one sonar beam is vertical, then the frequency spectra in the 
can be separated, and a crude estimate of direction obtained from the ratio of horizontal 
velocity spectra. But phase information is irrevocably lost through this procedure and 
the estimate is substantially biased when the waves are directionally spread. As a result, 
this estimator is not particularly useful, except perhaps in the case of swell. There is, 



however, phase information in the cross-correlations between the various range bins, 
and this fact allows the application of conventional signal processing techniques to 
estimate wave direction. 

The wave directional spectrum (WDS) is a mathematical representation of the 
5 wave direction as a function of azimuth angle and wave frequency, which is useful in 
describing the physical behavior of waves within the fluid medium. The most common 
existing devices used to obtain wave directional spectra are 1) pitch, and roll buoys, and 2) 
PUV triplets, described in further detail below. 

Pitch and roll buoys typically measure tilt in two directions as a surrogate for wave 

10 slope, along with the vertical component of acceleration. A recent variation uses GPS 
(Global Positioning System) measurements of three velocity components instead. The 
measured time series are Fourier transformed and the auto- and cross-spectra are formed, 
resulting in a cross-spectral matrix at each frequency. The elements of the cross-spectral 
matrix are directly related to the first five Fourier coefficients in direction (through 26) of 

15 the wave directional spectrum at each frequency (see Appendix Al). These buoys are 
typically used in deeper water. Unfortunately, the transfer functions for these buoys are 
complex, non-linear, and often difficult to determine. Additionally, the presence of a 
mooring line for the buoys adds additional complexity to the analysis due to added 
motion. Furthermore, such buoys are comparatively costly, vulnerable to weather and 

20 theft, and are not capable of measuring currents or wave heights. 

PUV triplets (so named due to their measurement of pressure and both 
components of horizontal velocity, namely u and v) are basically single point 
electromagnetic current meters having an integral pressure transducer. Time series of 
pressure and horizontal velocity from PUV triplets are processed in a manner similar to 

25 the measurements made by pitch and roll and GPS buoys, also giving only the first five 
Fourier coefficients in direction at each frequency. PUV triplets are typically bottom 
mounted, and generally only useful in shallow water. This significant disability is due to 
the decrease in high frequency response resulting from the decay of wave velocity and 
pressure with increased water depth. 
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Figure 1 illustrates a third and less common prior art technique for measuring 
wave directional spectrum employed by Krogstad, et al ( see "High Resolution Directional 
Wave Spectra from Horizontally Mounted Acoustic Doppler Current Meters," Journal of 
Atmospheric and Oceanic Technology, Vol. 5, No. 4, August 1988) as part of the 
5 CUMEX (Current Measurement Experiment) program. This technique utilizes an 
acoustic Doppler sonar system having a transducer array mounted on an underwater 
structure. The array is configured with sets of horizontally oriented acoustic transducers 
which project two acoustic beams in a horizontal plane 90 degrees apart. Beam 
propagation is therefore essentially parallel to the surface of the water, and skims the 

10 surface of the water as the beam disperses. Such a surface skimming geometry provides a 
relatively dense and uniform set of time lagged echoes and therefore permits estimation of 
the joint frequency-wavenumber spectrum S(f,k). See "Open Ocean Surface Wave 
Measurement Using Doppler Sonar " Pinkel, R. and J.A. Smith, J. Geophys., Res. 92, 
1987. Specifically, the directional spectrum D(Q) is expanded into a Fourier series, the 

15 coefficients of which are determined from the cross-spectral coefficient matrices generated 
from data obtained by the system. Since the acoustic beams are horizontal, no vector 
quantity (i.e., sensitivity vector) relating the beam geometry to the received current data is 
necessary. This technique is well suited to applications requiring only wave direction 
measurements and where a large, stable platform, such as a tower or large spar buoy, is 

20 available. However, there are a large number of applications, particularly in coastal 
oceanography and engineering, where it is desirable to know both the wave direction and 
vertical current profile, which the horizontal beam system can not provide. These 
applications include the analysis of sediment transport, atmosphere/sea interaction, 
pollutant dispersal, and hydrodynamic forces on off-shore structures. Additionally, it may 

25 be desirable in certain situations to simultaneously obtain wave height data along with the 
direction and current data. Due to the beam geometry, horizontal beam systems also are 
unable to measure current velocity above the wave troughs, which may be useful for 
studies of wave kinematics. 
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In summary, existing wave direction measurement techniques generally have 
several significant drawbacks (depending on type) including 1) inability to measure fluid 
current velocity and/or wave height along with WDS, 2) inability to readily measure wave 
directional spectrum at a broad range of depths; 3) inability to measure velocity profile 
5 above the wave troughs, 4) high degree of non-linearity; 5) high cost relative to the data 
obtained therefrom; and 6) susceptibility to damage/degradation from surface or near- 
surface influences. 

Accordingly, a system and method for accurately measuring the wave directional 
spectrum and current profile in a broad range of water depths is needed by both 

10 commercial entities and researchers. Such a system and method would further allow the 
optional measurement of wave height in conjunction with WDS. Additionally, the system 
would be highly linear in response, physically compact and largely self-contained in 
design, and could be deployed in a number of different scenarios such as being bottom 
mounted, moored, or even mounted on a mobile submerged platform. The flexibility of 

15 configurations permits the user to have the maximum degree of operational flexibility and 
device longevity without significantly impacting performance or accuracy. Additional 
benefits of economy would be obtained through the use of commercially available off-the- 
shelf broadband or narrowband sonar equipment where practical. 

20 Summary of the Invention 

The above-mentioned needs are satisfied by the present invention which includes a 
system and method for measuring the wave directional spectrum, current velocity within a 
given range cell or set of range cells, and wave height associated with a fluid medium by 
using acoustic signals. The present invention allows for accurate measurements of these 

25 parameters from fixed, moored, or mobile platforms using conventional Doppler sonar in 
conjunction with an upward and/or downward looking transducer array. 

In a first aspect of the invention, there is an improved system for measuring the 
wave directional spectrum of a fluid medium. In one embodiment, a broadband acoustic 
Doppler current profiler (ADCP) is used in conjunction with an upward looking, bottom 
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mounted multi-transducer (or phased) acoustic array to generate multiple acoustic beams, 
and sample a plurality of different range cells corresponding to different depths within 
those beams, in order to derive velocity and wave height data. Pressure readings are also 
optionally obtained from an associated pressure transducer located either on the ADCP or 
5 remotely thereto. This velocity, wave height, and pressure data is Fourier-transformed by 
one or more signal processors within the system (or in post-processing), and a surface 
height spectrum produced. A cross-spectral coefficient matrix at each observed frequency 
is also generated from this data. A sensitivity vector specifically related to the ADCP's 
specific array geometry is used in conjunction with maximum likelihood method (MLM), 

10 iterative maximum likelihood method (IMLM), iterative eigenvector (EEV), or other 
similar methods to solve the forward relation equation at each frequency and produce a 
wave directional spectrum. The ADCP may further be used to measure current velocity 
for various range cells in conjunction with the WDS and wave height measurements. 

In a second embodiment of the wave direction measuring system of the present 

1 5 invention, the sonar system and acoustic array are mounted to a platform, such as a vessel 
hull, with the array positioned so as to generate a plurality of upward and/or downward- 
looking acoustic beams such that altitude and bottom velocity, along with mean current 
profile, can be measured. 

In a second aspect of the invention, there is an improved algorithm and method of 

20 measuring the wave directional spectrum associated with a fluid medium. Specifically, a 
plurality of vertically oriented (possibly slanted) acoustic beams are generated using the 
system previously described. A plurality of different range cells within those beams are 
sampled to derive current velocity (and optionally, wave height) data. Pressure readings 
are also optionally obtained from the associated pressure transducer. Initially, the sampled 

25 data is processed to remove data outliers and calculate mean values for current velocity, 
depth, and pressure. Intrinsic wave frequency (f) and wavenumber magnitude (k) are also 
calculated for each observed frequency component. Next, a non-directional surface height 
spectrum is calculated by computing the power spectra for the range cells of interest, and 
the transfer functions for the sensor array. A cross-spectral correlation matrix is then 
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generated by selecting data from certain range cells, which may be the same or different 
from those selected for computing the power spectra above, applying a fast Fourier 
transform (FFT), cross multiplying all possible pairs of spectral coefficients for each 
observed frequency, and then averaging the results over time (from repeated Fourier 
5 transforms of sequential time segments) and/or over frequency (within bands of adjacent 
observed frequencies). After the cross-spectral matrix has been obtained, maximum 
likelihood, iterative maximum likelihood, and/or iterative eigenvector solution methods 
are applied to the array sensitivity vector, which is uniquely related to the chosen array 
geometry, and cross-spectral matrix to obtain a frequency specific wave directional 

10 spectra. Ultimately, the frequency-specific spectra are combined to construct a complete 
wave directional spectrum descriptive of both azimuth and frequency. The inventive 
algorithm can be run on existing signal processing equipment resident within the ADCP 
system, or run on an external computer as a post-processing task. 

In a third aspect of the invention, there is an improved method of measuring the 

15 wave height spectrum of a fluid medium using a submerged sonar system. A broadband 
ADCP sonar system of the type described above is used in conjunction with an upward- or 
downward-looking transducer array in order to measure wave height using one or a 
combination of methods. The first method is to determine the slant range to the surface 
using the measured backscatter intensity and/or signal correlation to interpolate the 

20 location of the surface. The wave height spectrum is determined as the power spectrum of 
the surface elevation. The second method is to measure beam velocity data from selected 
range cells within the beams and use the relationship between the velocity spectrum and 
the wave height spectrum from linear wave theory to determine the latter. Alternatively, 
in a third method, pressure measurements obtained from the ADCP's pressure transducer 

25 may be used to calculate wave height. 

These and other objects and features of the present invention will become more 
fully apparent from the following description and appended claims taken in conjunction 
with the accompanying drawings. 
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Brief Description of the Drawings 
Figure 1 illustrates a prior art wave directional spectrum measurement device 
utilizing two horizontal acoustic beams. 

Figure 2a is a perspective view of a first embodiment of an acoustic sonar system 
5 used for measuring wave directional spectrum, wave height, and current profile according 
to the present invention, the system being mounted to the bottom of the fluid medium 
volume. 

Figure 2b is a perspective view of a second embodiment of an acoustic sonar 
system used for measuring wave directional spectrum, wave height, and current profile 
10 according to the present invention, the system being mounted to the bottom via a mooring 
line or tether. 

Figure 2c is a perspective view of a third embodiment of an acoustic sonar system 
used for measuring wave directional spectrum, wave height, and current profile according 
to the present invention, the system being mounted on a submerged moving platform. 

15 Figure 2d is a perspective view of a fourth embodiment of a bottom mounted 

acoustic sonar system used for measuring wave directional spectrum, bottom velocity, and 
altitude according to the present invention, the system being mounted on a mobile surface 
platform in a downward-looking orientation. 

Figures 3 a and 3b are side elevational and top plan views, respectively, of an 

20 exemplary current profiler transducer array in the Janus configuration used for performing 
wave directional spectrum and wave height measurements in accordance with the present 
invention. 

Figure 3c is an overhead and perspective view of an exemplary transducer array in 
the "star" configuration optionally used in conjunction with the present invention. 
25 Figure 4 is a side view of the bottom mounted sonar system of Figure 2a 

illustrating the angular relationships of the acoustic beams to the surface of the fluid 
medium, and the relative position of range cells within the beams. 
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Figure 5 is a block diagram of an exemplary embodiment of the electronics for a 
broadband acoustic Doppler current profiler (ADCP), which may be used in an 
embodiment of the present invention. 

Figure 6 is a flow diagram illustrating the general method of determining wave 
5 directional spectrum according to one aspect of the present invention. 

Figure 6a is a flow diagram illustrating one exemplary embodiment of the initial 
data processing function of the method of Figure 6. 

Figure 6b is a flow diagram illustrating one exemplary embodiment of the wave 
height spectrum calculation function of the method of Figure 6. 
10 Figure 6c is a flow diagram illustrating one exemplary embodiment of the cross- 

spectral coefficient matrix calculation function of the method of Figure 6. 

Figure 6d is a flow diagram illustrating one exemplary embodiment of the 
maximum likelihood method (MLM) calculation function of the method of Figure 6. 

Figure 6e is a flow diagram illustrating one exemplary embodiment of the iterative 
1 5 maximum likelihood method (IMLM) calculation function of the method of Figure 6. 

Figure 6f is a flow diagram illustrating one exemplary embodiment of the function 
of constructing the two-dimensional wave directional spectrum D(Q,f) according to the 
method of Figure 6. 

Figure 7 is a plot of the linear wave dispersion relation obtained from an 
20 experimental deployment of the WDS measurement system using an embodiment of the 
present invention. 

Figure 8 is a mesh plot of an exemplary maximum likelihood method(MLM) 
estimate of the directional wave spectrum generated using an embodiment of the present 
invention. 

25 Figure 9 is a wave height spectrum obtained by integrating the spectra of Figure 8 

over all azimuthal angles. 

Figure 10 is a plot comparing estimates of wave directional spectrum versus 
azimuth for various input spectra as derived using both the maximum likelihood and 
iterative maximum likelihood methods. 
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Detailed Description of the Preferred Embodiments 
Reference is now made to the drawings wherein like numerals refer to like parts 
throughout. 

5 Figure 1 shows a transducer array for measuring wave height as is known in the 

technology. Figure 2a illustrates a first exemplary embodiment of the WDS measurement 
system of the present invention. As shown in Figure 2a, the system 100 is bottom- 
mounted and includes a body element 101 containing sonar electronics and processing 
equipment 102, and a multi-transducer array 103 having the individual transducer 

10 elements arranged in the Janus configuration (see discussion of Figures 3a and 3b below). 
This transducer array 103 generates acoustic beams 104 which are coplanar in the vertical 
plane 106 yet divergent from a horizontal plane 108 parallel to the surface of the fluid 
medium 1 10. The fluid medium 1 10 is most often natural or man-made bodies of water, 
especially the ocean. It should be noted that while the Janus array configuration is used in 

1 5 the embodiment of Figure 2a, other array configurations which form beams having an 
angular relationship to the horizontal plane 108 may also be used. For example, a 
"pinwheel" array (e.g., one where the acoustic beams are skew-divergent from the 
longitudinal axis of the array), or "star" array (non-coplanar, non-skewed beams) may also 
be used. Additionally, phased or time-delayed arrays may be used in conjunction with the 

20 present invention. Note that the array sensitivity vector H, which is described in greater 
detail below, is unique to the specific array geometry employed, and accordingly is 
modified for use with such alternate array configurations. 

Referring again to Figure 2a, the sonar system 100, including the transducer array 
103, is mounted to the bottom 110 using any number of well known bottom mounting 

25 techniques including attachment to a heavy frame, trawl-resistant cage, or buried post. 
Alternatively, the system 100 can simply be laid on the bottom (or partially buried) in the 
desired orientation with no mounting hardware or attachment; a significant improvement 
over prior art systems which require a mooring line or platform. Since there is effectively 
no motion of the array within local spatial frame of reference 1 12 when the system 100 is 
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bottom mounted, no errors resulting from the movement or rotation of the array 103 are 
generated, and hence the accuracy of the system is maximized. Additionally, bottom 
mounting places the system 100 away from the potentially damaging influences of surface 
craft, strong surface waves and current, and vandals. For these reasons, fixed mounting to 
5 the bottom 1 10 as depicted in Figure 2a is the preferred method; however, the system 100 
may alternatively be mounted via a mooring line (Figure 2b), or on a moving underwater 
platform (Figure 2c), and still produce useful WDS, current, and wave height data. 

As shown in Figure 2d, the transducer array 103 may also be inverted for use with 
a surface or near-surface application, such as within the hull of a surface vessel 120, such 

10 that downward-projecting acoustic beams 104 are generated. In this way, WDS or current 
velocity at varying depths and at the bottom 1 10 of the fluid volume can be measured. 
Array altitude above the bottom (e.g., the height of the array above the local bottom), 
which may be different than water depth, can also be measured using this configuration. 

Referring now to Figures 3a and 3b, an exemplary transducer array configuration 

15 is shown. As previously discussed, the Janus configuration is preferred for use with the 
present invention, although other transducer configurations may be used. The Janus 
configuration utilizes transducer elements 140 which are spaced at 90 degree intervals in 
the azimuthal plane 142, and which are angled with respect to the longitudinal axis 144 of 
the array. The construction and operation of the Janus configuration array are well known 

20 in the sonar technology, and accordingly shall not be further described. 

Figure 3c illustrates a "star" (aka "bugeye") transducer configuration which can be 
used as an alternative to the Janus configuration previously described. As shown in the 
Figure, the star array generates four separate acoustic beams: three beams inclined with 
respect to the longitudinal axis 144 of the array and set at 120-degree intervals in azimuth, 

25 and a fourth beam centrally located and coincident with the longitudinal axis 144. 

As is discussed in greater detail with respect to Figure 5 below, an acoustic 
Doppler current profiler (ADCP) is the preferred sonar system used in the wave 
directional measurement system 100 of Figures 2a-2d. Unlike prior art WDS 
measurement systems which must be located within or close to the strata of the fluid 
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medium being analyzed, the ADCP, using its upward-looking array, is capable of 
remotely sensing current velocity and wave height, thereby providing the WDS 
measurement system 100 with greater range capability as well as protection from surface 
and near-surface hazards. 
5 Figure 4 is a side view of the bottom mounted sonar system 100 of Figure 2a, 

including the transducer array of Figures 3a and 3b, illustrating the angular relationships 
of the acoustic beams 104 to the surface of the fluid medium, and the relative position of 
individual range cells within the beams. As shown in the Figure, the beams project 
upward from the transducer array 103 at an angle, typically between 20 and 30 degrees, 

1 0 although other values may be used, in relation to the vertical axis 200 of the local frame of 
reference 112. It is assumed for illustration that the longitudinal axis of the array 144 is 
coincident with the vertical axis 200 of the local frame of reference 1 12. However, it can 
be appreciated that the array 103 may be mounted at somewhat of an angle ("tilt") with 
respect to the vertical axis 200 depending on bottom type, contour, and other factors. 

15 Such a tilt is accounted for by selecting different range cells or bins within the various 
beams, and using the precise location of each range cell in the calculation of the 
sensistivity vector. Additionally, when the array 103 is mounted on a mooring line (as in 
Figure 2b) or on a submerged moving platform (Figure 2c), rotations off of the vertical 
axis 112 may be induced. Depending on the magnitude of the rotation, effects on the 

20 ultimate calculation of WDS, wave height, and current velocity within a given set of range 
cells may exist. Accordingly, a correction algorithm may be applied if desired to 
compensate for such movement. 

Figure 5 illustrates an exemplary embodiment of the electronics for a broadband 
acoustic Doppler current profiler (ADCP) such as a Rowe-Deines Instruments Model 

25 BBADCP VM-150 used within the present invention. While the following discussion 
may refer to this ADCP system, it can be appreciated that other models and types of sonar 
systems, such as narrowband Doppler systems or non-Doppler-based systems, may be 
used with the present invention depending on the particular application and needs of the 
user. 
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Referring again to Figure 5, the transducer array 103 is electrically connected to 
the electronics assembly 170 which includes a mixer network 172, low pass filter network 
174, sampling module 176, and digital signal processor (DSP) 178. Signals generated by 
the transducer array elements 140 upon the receipt of acoustic signals are fed via the 
5 transmit/receive switches 180 to preamplifiers 182 and receiver amplifiers 184 which 
condition and amplify the signal(s) for further processing by the electronics assembly 170. 
A coder transmitter 186 and power amplifier 188 are used in conjunction with the DSP 
178 to feed transmission signals to the transducer elements 140 via the transmit/receive 
switches 180. Thus, the same transducer elements are used for both transmit and receive 
10 functions. Additional details regarding the exemplary broadband ADCP system are 
contained in U.S. Patent No. 5,208,785, "Broadband Acoustic Doppler Current Profiler" 
assigned to Rowe-Deines Instruments, Inc., which is incorporated herein by reference in 
its entirety. 

It should be emphasized that the above-described system 100 utilizes a standard 
15 Rowe-Deines or other comparable ADCP, and requires no special hardware or 
modifications. This fact greatly increases the design and manufacturing economy of the 
present invention, and allows for easy retrofitting to existing ADCP systems. A pressure 
sensor is used to measure mean water depth (as well as optionally the wave height 
spectrum as discussed further below. However, such sensors are often available as options 
20 on commercial ADCP systems, or can be otherwise easily integrated into or mounted on 
the ADCP. Most any commercially available sensor adapted for underwater use, such as 
the Sensym Hastelloy C22 flush mount pressure sensor cell, can be used in this 
application. 

25 Method and Algorithm for Estimating Wave Directional Spectrum and Wave Height 

Each of the previously described hardware embodiments of the wave measurement 
system of the present invention utilize a specially designed algorithm to calculate WDS 
and wave height. This algorithm is ideally embodied in the form of software running on 
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the existing signal processing capability of the ADCP (or other chosen sonar system), 
although it can be appreciated that the algorithm or components thereof could be 
embodied in firmware and/or hardware as well. The algorithm is based on a general 
method of calculating WDS which uses "sensor" data (e.g., current velocity, wave height, 
5 and/or pressure data) in conjunction with a sensitivity vector H, the latter being uniquely 
related to the chosen array geometry. The following discussion provides a theoretical 
overview of this method, followed by a detailed description of the operation of the 
algorithm used in conjunction with the preferred ADCP system described above. 
Appendix I provides additional detail on the derivation of the sensitivity vector H, cross- 
1 0 spectral matrix C, and wave directional spectrum D(Q , /) . 

Theoretical Overview 

ADCP sonar systems measure, inter alia, the instantaneous current velocity 
component projected along the axis of each of its angled acoustic beams. The current 

15 can then be computed as the mean of the difference in velocities between opposing 
beams. Since successive positions along the angled beams correspond to different 
horizontal locations, the set of range cells within the beams constitutes a spatial array. 
Useful information regarding wave direction is contained in the current velocity cross- 
spectra (e.g., the array covariance matrix). Several other factors relating to the acoustic 

20 beam geometry must be considered, however. 

First, the velocity signal-to-noise ratio (SNR) varies as a function of depth. At a 
given frequency, this change is due largely to the decay of wave velocity and pressure 
with depth. Hence, for the shorter waves, only range cells fairly close to the are of 
practical use. 

25 Second, velocity measured by the sonar is a linear combination of both vertical 

and horizontal wave velocities, the relative weighting therebetween being a function of 
both wave direction and water depth. A mathematical relationship connecting the 
along-beam component of wave velocity and surface elevation is needed. 
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For practical purposes, a random wave field can be treated as if it were a linear 
superposition of infinitesimal uniform plane waves, each with a particular directional 
orientation, frequency, and wavenumber. Neglecting non-linear effects, the dispersion 
relation is used to constrain these plane waves to a set that satisfies the linear dispersion 
5 relation connecting frequency and wavenumber. When there is no mean current, the 
linear dispersion relation for an infinitesimal plane wave is: 

f = — ^jgktanh(kh) Eqn. 1 

2k 

where /is the wave frequency (Hertz), 

g is the gravitational acceleration constant (about 9.8 m/s , 
10 h is the water depth (meters), 

k - 2nfL is the wavenumber magnitude (radians per meter), and 
L is the wavelength (meters), 
Allowing for a steady uniform current u c and assuming an observation platform fixed to 
the bottom, the linear dispersion relation has the somewhat modified form: 

.5 / -' / + ^"-- k Eqn. 2 

= £ [ ^gk tanh(Jtt) - u c k cos(p -9) ] 

where f obs is the observed wave frequency (Hz), 

/ is the intrinsic wave frequency that would be observed in a reference frame 

moving with the current (from Eqn. 1), 
k is the wavenumber vector perpendicular to the wave crests (radians/m), 
20 pointed in the direction the wave moves relative to the water, 

P is the azimuth angle of the direction the current is moving toward, 
G is the azimuth angle opposite to the direction of the wavenumber vector k, (the 

direction the waves are coming from relative to the water), 
u c is the current vector, and 
25 u c = |u c | is the current speed (m/s). 
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Note that velocity contributions of order (ak) 2 and higher are neglected, where a is a 



measure of wave amplitude (e.g. ak is proportional to wave slope). Terms of order (ak) 2 
in the dispersion relation are also discarded. Since the waves of interest are those near 
the spectral peak, which conservatively have slopes less than 0.1, these corrections are 



If the current is not uniform with elevation, the appropriate value of u c to use at 
each intrinsic frequency is the weighted-average velocity profile, using a normalized 
weighting function proportional to the square of the horizontal wave orbital velocity: 



where <u(z)> is the time-averaged horizontal current profile. 

For the wave directional spectrum, the usual convention is to use the linear 
dispersion relation to eliminate the magnitude k of the wavenumber as an independent 
15 variable. The wave directional spectrum D(Q,f) therefore represents the power 
spectral density of the infinitesimal plane waves in the two-dimensional azimuth- 
frequency space. 

In general, an instrument for measuring the wave directional spectrum can make 
use of a combination of velocity, pressure, and/or surface elevation measurements at a 

20 single point ("triplet") or at an array of points at or below the water surface. Linear 
wave theory can be used to relate the directional wave spectrum D(Q,f) to these 
measurements. The Fourier transforms of measurement time series can be cross- 
multiplied to form cross spectral coefficients, which can be arranged in a cross-spectral 
coefficient matrix C at each observed frequency. The so-called "forward relation" is a 

25 theoretical model connecting the unknown wave direction spectrum D(Q 9 f) with the 
observable array covariance matrix C(f obs ) (the matrix of velocity cross-spectra between 
the various sonar range cells). The forward relation takes the following form: 



5 



less than 1%. 



10 




Eqn. 3 
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C(f 0 J= f m,f oks )D(Q,f)HHQ,f obs )JdQ 



Eqn. 4 



where f obs =f obs (f) is the mapping given by Eqn. 2, 

H is the sensitivity vector (known from linear wave theory), 
H t is the complex conjugate transpose (Hermitian) of H, 




9 



1- 



u c cos(p-e) 



c 



8 



is the Jacobian, 



Eqn. 5 



d(2nf) 
dk 




sinh(2£/*) 



2kh 



is the group velocity, and 



Eqn. 6 



c 



P h 




— tznh(kh) is the phase velocity 



Eqn. 7 



Note that the form of the sensitivity vector H is related to the geometry of the array 
formed by range cells selected from those in the ADCP beams. Each selected range cell 
from a particular beam and depth corresponds to an element of H. The sensitivity 
vector H may also optionally contain elements corresponding to surface elevation 
measured acoustically at each beam, and/or a pressure measurement. 

Note also that when there is a current, there may be more than one / that maps to 
eachy^ 5 , in which case C is the sum of such mappings. Since negating the frequency is 
equivalent to reversing the direction o the wavenumber, 6 should be replaced with 0+tt 
whenf obs is negative due to strong current. 

Again, the forward relation (Eqn. 4) applies to any array or triplet that measures 
wave directional spectra. The inversion of the forward relation to determine wave 
directional spectrum from measured cross-spectral matrix C presents significant 
complexities in that an infinite-dimensional object (D) is to be estimated using a set of 
finite-dimensional observations (C). Hence, the problem is severely underconstrained and 
no unique solution exists. A number of different solution techniques have been developed 
to address the inversion problem. Of most use to the present invention are the so-called 
maximum likelihood method (MLM) and iterative maximum likelihood method (IMLM), 
the general application and operation of which are well known in the signal processing 
arts. See Appendix I for a more detailed discussion of the theoretical foundation for these 
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techniques. It should be noted that while the present invention is described in terms of the 
MLM and EMLM techniques, other techniques such as the iterative eigenvector (IEV) 
method, may be applied depending on the particular application. 



5 Description of WDS Algorithm 

Referring now to Figure 6, the algorithm for determining the wave directional 
spectrum according to the present invention is described. In the present embodiment, the 
algorithm 300 consists of the following general steps: 1) conducting initial data processing 
302; 2) calculating a non-directional wave height spectrum 304; 3) calculating a cross- 

10 spectral matrix C for each observed frequency 306; 4) estimating the wave directional 
spectrum (for each observed frequency) 308a, 308b; and 5) constructing a complete two- 
dimensional wave directional spectrum from the estimates derived in Step 4) 309. Each 
of these steps is described in detail below, with reference to Figures 6a-6f, respectively. It 
should be noted that the individual steps 302, 304, 306, 308a, 308b, 309 of the algorithm 

15 300 (and their sub-steps) may be permuted under certain circumstances. Furthermore, 
while specific computations are described, it can be appreciated that other mathematical 
approaches or techniques may be substituted. 

As shown in Figure 6a, the initial processing step 302 of the algorithm 300 
includes a series of sub-steps 310, 312, 314, 316, 318, 320 which collectively prepare the 

20 raw data obtained from the sensors for further processing. In the first sub-step 310 of 
algorithm step 302, raw data obtained from the sensors is decoded using any number of 
decoding algorithms known in the prior art. In the second sub-step 312, the surface height 
for each acoustic beam is calculated. This is accomplished using the interpolated peak 
location in the acoustic backscatter intensity, as described in "Measuring Wave Height and 

25 Direction Using Upward-Looking ADCPs", Terray, E., et al, IEEE Oceans 1997, August, 
1997, which is incorporated herein by reference in its entirety, although other methods 
may be used. The third sub-step 314 involves processing the decoded raw data to identify 
"outliers," data having a poor correlation, or data associated with other anomalies present 
in the fluid medium (such as marine life). Next, in the fourth sub-step 316, "bad" data 
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identified in sub-step 314 is replaced with data derived by interpolating between "good" 
data values. Mean values of current velocity (and optionally water depth and pressure) are 
then calculated in sub-step 318 from the processed data resulting from the data 
replacement process 316. Finally, in the last sub-step 320, values of wave frequency (f) 
5 and wave number magnitude (k) are calculated using numerical inversion of the linear 
dispersion relation (see eqn. 2 above, which is discussed in additional detail in Appendix 
I) 

As shown in Figure 6b, the step of calculating the wave height spectrum S H 304 in 
Figure 6 includes a series of sub-steps 322, 324, 326, 328, 330, 332, 334 which 

10 collectively calculate S H from the data initially processed in the first algorithm step 302. 
In the first sub-step 322 of algorithm step 304, velocity data is selected from certain range 
cells within the acoustic beams 104. Using the apparatus of Figure 2a above, these range 
cells are selected primarily on their proximity to the surface of the fluid medium (so as to 
reduce the effects of wave velocity and pressure decay with depth, as previously 

15 described), although other cells may be chosen. Optionally, surface height and pressure 
data is also selected. It should be emphasized that while the present algorithm 300 has the 
capability to incorporate wave height and pressure measurements within the sensitivity 
vector H (see discussion of Figure 6d below), such measurement are not required to 
calculate the WDS. The use of such measurements do, however, generally increase the 

20 relative accuracy of the estimate provided by the system. 

In the second sub-step 324 of Figure 6b, the selected velocity, surface height, and 
pressure data is detrended by subtraction of a trendline determined by least-squares fit and 
windowed by multiplication by a window function oi time such as a Bartlett window, for 
the purpose of reducing spectral leakage. The third sub-step 326 computes a power 

25 spectrum for each set of data by performing a fast Fourier transform (FFT) and squaring 
the magnitude of the result, although other methods may be used. Next, in the fourth sub- 
step 328, the power spectra are averaged into frequency bins (so-called "decimation"). 
Transfer functions T(f obs ) = |Hf 2 are then computed for each "sensor" in sub-step 330 
using Eqns. 16 and 17 below. It should be noted that the term "sensor" in the present 
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context refers to either 1) current velocities obtained from the range cells of the acoustic 
beams at the same water depth; 2) the surface height measurement obtained for each 
acoustic beam (see calculation of surface height sub-step 312 above), or 3) the pressure 
sensor. Combinations of these sensors may also be used to calculate the wave height 
5 spectrum. 

In the sixth sub-step 332 of Figure 6b, power spectra for groups of range cells at 
the same depth are combined by adding the power spectral coefficients. Finally, in the last 
sub-step 334, all power spectra previously calculated are combined into one wave height 
spectrum S H using, for example, a least squares fit technique and the transfer functions T(f 

10 obs). 

Referring now to Figure 6c, the method of calculating of the cross-spectral matrix 
is described in detail As in the calculation of S H described above (Figure 6b), data from 
specific range cells, surface height data, or pressure sensor data are selected in the first 
sub-step 336. Note, however, that the data selected in this sub-step 336 may be the same 

15 or different from that selected as part of the calculation of S H 322. Smaller range cells may 
be needed to resolve the direction of the shorter waves, while range cells may be 
combined into larger ones (by averaging velocity measurements) to reduce the noise level 
for computing the wave height spectrum. 

Again, the selected data is detrended and windowed (sub-step 338). After 

20 windowing and detrending, complex spectra are computed (using FFT or another 
comparable technique) 340. Next, all possible pairs of spectral coefficients for each 
observed frequency (f obs ) are cross-multiplied 342. In order to perform this cross- 
multiplication sub-step 342, one spectral coefficient from each pair is first complex- 
conjugated. Lastly, the cross-spectral matrices for each observed frequency are averaged 

25 over time (and/or within frequency bands) to produce the cross-spectral matrix as a 
function of observed frequency 343. 

Referring now to Figure 6d, the maximum likelihood method (MLM) step 308a of 
Figure 6 is described in detail. As previously described, the MLM is one algorithm 
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available for solving the so-called "forward relation" which relates the wave directional 
spectrum D (Qf) with the observable cross-spectral matrix C (f obs ). 

As shown in Figure 6d, the first sub-step 344 of the MLM is to invert the cross- 
spectral coefficient matrix C for each observed frequency using a singular value 
5 decomposition (SVD) technique. SVD techniques are well known to those skilled in the 
signal processing arts, and accordingly will not be discussed in further detail herein. 

Next, the array-specific sensitivity vector H is calculated 346. Each selected range 
cell from a particular acoustic beam and depth corresponds to an element of H. Note that 
H may also contain elements corresponding to optional surface height and pressure 
10 measurements. The sensitivity vector H represents the ideal measured response to a 
plane wave, assuming linear theory and the linear dispersion relation (see previous 
discussion), and neglecting instrument noise. The linear wave model predicts: 



T| (jc, /) = a cos(k • (x - u c t) - 2k ft) -a cos(k • x - 27r f obs t) 



Eqn. 8 



15 



p(x,z,t) = pga 



cosi 



ih[k(h + z)] 



cos(k*x-2n f obs t)-pgz 



Eqn. 9 



cosh(M) 



u(x,z,t) = a(2n f) 



cosh[A:(A + z)] 
sinh(£/*) 




Eqn. 10 



w(x,z,/) =a{2n f) 



sinh[A;(/* + z)] 
sinh(^) 



sin(k«x-27i/ o6j 0 



Eqn. 11 



20 



where r| is the surface elevation (m), 
p is the pressure (Pascal), 
p is the water density (kg/m ), 
g is the gravitational acceleration (m/s 2 ), 
u is the horizontal velocity vector (m/s), 
w is the vertical velocity component (m/s), 
a is the half-amplitude of the plane wave (m), 



25 
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h is the water depth (m), 

x is the horizontal position vector, and 

z is the vertical position {-h at the bottom, o at the surface) (m). 
The Fourier transforms in time of the above sine waves are delta function line 
spectra, the areas under which are the respective Fourier coefficients (indicated by 
tildas): 

tT( x > fobs ) = l a ex P(* • x ) Eq n - 12 



10 



cosh[A;(/2 + z)] 
cosh(£/z) 



exp(/k • x) 



u(x, z, f obs ) = £ a(27r f) — .\\ 9 „ exp(*k • x) 



smh(kh) 



k) 



Eqn. 13 



Eqn. 14 



*, /«. ) = -| H2n f) + * )] exp(/k * x) Eqn. 15 

sinh(£A) 

Note that the Fourier coefficient of the surface displacement at x = 0 is simply a/2, 
15 regardless of the direction or frequency of the plane wave. The squared magnitude of 

this coefficient is the non-directional power spectrum of the wave height. Therefore the 

Fourier coefficients of the various measurements are scaled by 21a to represent the 

response to a plane wave of unit spectral density. 

For velocity measurements, each element of H is computed as the scaled Fourier 
20 coefficient corresponding to component of the velocity in the direction of the beam at 

the location of the range cell: 



H(z,x„) = - 



U + W1. 



•b. 



27r/exp[/k«x fl (z)] 



sinh(&/0 



cosh[k(h + z)] — ^ -isinh[k(h + z)] (i 2 »b„ ) 



Eqn. 16 
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where b n is the unit vector pointing outward in the direction of the n ih beam of the 
sonar system, 

x n (z) is the horizontal displacement vector for the range cell (m), and 
i z is the unit vector in the vertical direction. 



Note that the x origin should be chosen at the sonar system 100 (ADCP in the preferred 
embodiment) so that the horizontal component of the beam unit vector b n points in the 
same direction as x n . The mapping among k, f, and f obs are given by the linear 
dispersion relations and their inverses, which are calculated numerically. 
10 For optional surface elevation measurements at each beam, the corresponding 

elements of H are computed as the phase adjustments needed to account for the effect of 
the horizontal displacements of the beams from the x origin: 

2 

H(x n ) = - n = exp(zk • x„ ) Eqn. 1 7 

a 

where b n is the unit vector pointing outward in the direction of the n h beam of 
15 the ADCP. 

For the optional pressure measurement at the ADCP, the corresponding elements 
of H are: 

ur* ^ 2 ~ ^ coshfc^ + z^)] 

H(z ADCP ) = -p = pg — — Eqn. 18 

a cosh(Art) 

Once the sensitivity vector H is calculated per sub-step 346, the non-normalized wave 

20 directional spectrum D (0) is calculated per sub-step 348 (for the specific observed 

frequency under consideration) using the following MLM estimator: D(Q ) 

D(Q)= !* D t Eqn. 19 

H + C l U 

where H f is the complex conjugate transpose of the sensitivity vector H, and C" 1 is the 
25 inverse of C. In the final sub-step 350 of Figure 6d, the wave directional spectrum D is 
normalized such that: 



-22- 



2n 



jDdQ =1. 

0 

Referring now to Figure 6e> the iterative maximum likelihood method (IMLM) 
step 308b of the algorithm of Figure 6 is described. Note that this step 308b is optional, 
and may be bypassed for one or more observed frequencies as desired. As shown in 
5 Figure 6e, the IMLM step 308b begins with the first sub-step 352 which is initialization 
with the normalized wave directional spectrum D n derived for each observed frequency in 
prior step 308a. (Note that the subscript "n" on the wave directional spectrum D is an 
index used to denote the iteration number. Hence, when n = 0, D n = D MLM ). Next, the 
forward relation (eqn. 3 above) is used to compute the iterative cross-spectral coefficient 
10 matrix C n using the following relationship (sub-step 354). 

C n = |HZ) n HVe Eqn. 20 

0 

The MLM step 308 of the algorithm 300 is again repeated in sub-step 356 using the 
computed value of C n in order to derive the matrix M, as follows: 

Af(£>.)= . K t Eqn. 21 

n> H t C -l H 

1 5 where K is chosen such that: 



2n 



jM(D n )dQ =1 Eqn. 22 

o 

Next, in sub-step 358, the index n is incremented by one (i.e., n+1) and the directional 
spectrum calculated according to the relationship: 

A, +1 = D n -co „ [M(D n )-D MLM ] Eqn. 23 

20 where o> R is the relaxation parameter, which is chosen to be approximately 1.1 in the 
present embodiment. 

In the next sub-step 360, the spectrum is fixed by zeroing negative values (or a 
similar procedure) to ensure that D n+1 is greater than or equal to zero. Subsequently, the 
spectrum is renormalized in substep 362 so that: 
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\*D n+] dQ=l Eqn.24 

Lastly, the convergence criteria is computed, (substep 364). Note that in the present 
embodiment of the algorithm 300, the convergence criteria are applied to each observed 
frequency (e.g., upon satisfying the IMLM convergence criteria for a given observed 
5 frequency, the MLM step 308a and optional IMLM step 308b are again performed for the 
next observed frequency). An exemplary set of criteria are as follows: 

i. a fixed maximum number of iterations (typically 3 to 5) has occurred 
without meeting any other stopping criteria; 

ii. the mean absolute difference in the smeared spectra integral 

1 0 | |(M (D n ) - D MLM ) dQ is greater than it was on the previous iteration; 

iii. the relative squared difference in the measured and forward-relation 
cross-spectral matrices is less than a predetermined threshold value: 



(c„-c„) 2 



< threshold value Eqn. 25 



where a is the estimated standard deviation of the i 9 j - th element of the measured C; 
15 or 

iv. the relative difference in the measured and predicted cross-spectral 
matrix C (Eqn. 25) is greater than it was on the previous iteration. 

The above-listed convergence criteria are applied such that satisfaction of any one 
(or more) of the criteria will terminate further iteration (for the observed frequency under 
20 analysis). When all observed frequencies have been analyzed, the algorithm then proceeds 
to the final algorithm step 309. 

As previously noted, the convergence criteria described herein are but one 
embodiment, and other criteria or mechanisms for terminating the IMLM iteration process 
may be substituted depending on the specific application of the algorithm and the needs of 
25 the user. 
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The final step 309 in the algorithm of Figure 6 is now described. As shown in 
Figure 6f, the complete two-dimensional wave directional spectrum is constructed by 1) 
renormalizing the observed frequency-dependent wave directional spectrum 370; 2) 
determining the presence of a non-zero current 372, and 3) in the event of non-zero 
5 current, remapping the frequency and applying the Jacobian operator 374. 

Renormalization of Z)(0 , f obs ) is accomplished by multiplying by S H so that: 

\D®,f obs )dQ=S H (J obs ). Eqn.26 

0 

Next, the magnitude of the current is examined, and if no current is present (magnitude = 
0), no further analysis or computation is performed, and the calculation of D(Q,f ob$ ) is 
10 complete. If, however, a non-zero current is present, the frequency is remapped using 
Eqn. 2, either by interpolating to a regular grid or plotting with a distorted frequency axis 
(that is, non-linear in observed frequency, linear in intrinsic frequency), and the Jacobian J 
is applied such that: 

15 D(Q,f) = D(Q,f obs {f})J Eqn. 27 



Where 



J 



obs 



df 



Eqn. 28 



20 Referring now to Figures 7 through 10, output from the aforementioned algorithm 

300 at various stages of the calculation process is illustrated. 

Figure 7 is an exemplary plot of the linear wave dispersion relation (eqn. 2) 
obtained from an experimental deployment of the WDS measurement system of the 
present invention. 

25 Figure 8 is a mesh plot of an exemplary maximum likelihood method (MLM) 

estimate of the directional wave spectrum generated using the present invention. 
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Figure 9 is a wave height spectrum obtained by integrating the spectra of Figure 8 
over all azimuthal angles. 

Figure 10 is a comparison of estimates of wave directional spectrum versus 
azimuth for various input spectra as derived using both the maximum likelihood and 
5 iterative maximum likelihood methods, using a known directional spectrum to simulate a 
cross-spectral matrix C. Note the improved performance of the IMLM estimate as 
compared to the MLM estimate. 

Calculation of Wave Height 

10 

As previously described, the wave directional measurement system of the present 
invention is further capable of calculating the significant wave height asssociated with 
waves in a the fluid medium. Specifically, the significant wave height Hs (measured in 
meters) is calculated as 4 times the square root of the area under the wave height spectrum 
1 5 SH (see previous discussion of Figure 6b for derivation of wave height spectrum) over the 
frequency range of interest, typically 0.03 to 0.5 Hz. This frequency range is expected to 
be dominated by waves rather than tidal fluctuations or noise. This calculation is similarly 
performed using the processors of the ADCP and/or as a post-processing task on an 
external computer. 

20 While the above detailed description has shown, described and pointed out the 

fundamental novel features of the invention as applied to various embodiments, it will be 
understood that various omissions and substitutions and changes in the form and details of 
the device illustrated may be made by those skilled in the art, without departing from the 
concepts of the invention. It should be further noted that while the present invention has 

25 been described in terms of a device and method useful in an underwater environment such 
as an ocean, harbor, or lake, it may also find application in other types of fluid mediums. 
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* 

Appendix I 



Al. The Array Response to Waves 

The sonar system 100 of Figure 2a is presumed to be oriented as shown in 
5 Figure 4, with the beams sequentially numbered from 1 to n, (n = 4 in Figure 2a) and 
having a coordinate system with z positive upwards (toward the surface, so that the fluid 
medium occupies the region 0 > z > -d. The x-axis is assumed to lie along beams 1 and 
3, with positive x pointing outwards along beam 1. Beams 2 and 3 respectively lie 
along the positive and negative y-axis. The coordinates of a range cell satisfy 

10 z = yjx 2 +y 2 I tana , where a is the beam angle with respect to the vertical 

We denote the component of velocity along the X-th beam (k = l,...n) by V^(x,0 
(where x is taken to lie at the center of a range bin); u{xj) denotes the physical wave 
velocity vector at location x. Taking a Fourier transform of these variables with respect 
to time yields the amplitudes V(x,<o) respectively. 

15 The relation between these amplitudes can be written as 

^ife<») 
F 2 (x,co) 
F 3 (x,(fl) 

y 4 (x 9 (o)_ 

or more compactly in matrix form as 

V(x y ay) = A(a)u(x y (o) A1.2 

The 4x4 cross-spectral matrix 
20 C(x, y | © ) =< F(x,co )V * (x> ) > A 1 . 3 

must ultimately be related to the wave direction spectrum D(co,0) . Note that V(x 9 to) 
is the Fourier amplitude at angular frequency co> "t" signifies the complex conjugate 

transpose, and <...> denotes a statistical average over realizations of V . 
Substituting equation A1.2 into A1.3 gives 



No menclature may vary from that used in "Detailed Description" section 

-27- 



sina 
sin a 
0 
0 



0 cosa 

0 cosa 

sina cosa 

-sin a cosa 



M 2 (X,C0) 

« 3 (x,cd)_ 



AL 1 



C(x y x^_ | co ) = A < u (x , co )u + (x^, co ) > A ' A 1 . 4 

where "f" denotes the transpose. 

For linear waves propagating in water of depth d, u is given by 

«(*,©) = Jrf 2 *II(*,a) Al. 5 

5 where r}(£,co) is the Fourier amplitude of the surface displacement field at angular 
frequency co and wavenumber k ~ (k x ,ky), 

m,®\x) = je ik -Hk x F,k y FHkF'l M g 

F(kz) = cosh(k(z + d)) I sinh(JW), and F'(x) = dF/dx. 

Then C becomes 

C(x,x'|co) = \d 2 kA\W{k^)ll-A^ A1.7 

1 0 where ^(^co) denotes the frequency/wavenumber spectrum of r\ (x, t) . 

Equation A1.7 is the forward relation between C and and defines a 3- 
dimensional inversion problem. This can be reduced to a 2-dimensional estimation 
problem by using the dispersion relation, defined by 

co 2 =g*tanh(fc/), k=\k\ Al. 8 

1 5 where g = 9.8 m/s is the acceleration due to gravity. 

This relation can be used to define frequency and frequency/direction spectra 
(denoted by S(<$) and £>(co ,9 ) respectively) via 

co )d 2 k = 5 (co - a )S(o )D(a ,9 )dcdQ Al . 9 

where 8 is the Dirac delta function, and a(k) satisfies the dispersion relation A1.8 for 
20 each wavenumber L The direction spectrum is normalized to unity area 

I d9D(co,9) = l ALIO 

Substituting into equation A1.7, we have 

C(x 9 x[\(o)= rddH^ 9 <Q\^)D(to 9 Q)H f (^ t o)\x) Al. 11 

where H(Q ,© | jc) = A(a)U(k(®) 9 (o \ x)> and k(<$) satisfies the dispersion relation for 
25 frequency co. 
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Since the wave field is stationary then (except for some short-range correlations 
introduced by windowing) the various spectral bands in frequency are uncorrelated, and 
equation ALU can be considered to hold for each frequency separately. 

Finally, for later purposes, it is easier to consider equation Al.ll as a single 
5 matrix equation at each frequency, rather than as a 4x4 matrix equation for each spatial 
separation x-tf. So we assume that the N range cells along each of the four beams 
have been ordered in a list (the order is unimportant). Then C becomes a 4N x 4N 
matrix, and H becomes 4N x 1. Note that C is still complex conjugate symmetric, and 
so only 2N(4N+1) elements are independent. 
10 When ordered in this way, the dependence on x and x!_ can be omitted, and 

equation Al.ll written as 

C= l + *ddH(Q)D(Q)H\8) Al. 12 

where it is understood that a separate equation for each frequency is implied. Note, if 
we consider a plane wave propagating in the direction 0 O , then Z>(0) = 8(9 -9 J , and 

15 C = H(Q 0 )H 1 (9 0 ) , so that H gives the directional response of the array. 

Equation A1.12 connects a finite set of data (contained in C) to the continuous 
function D(Q ) . This expression can be converted into a matrix equation suitable for 
parameter estimation by approximating D by a finite set of coefficients. As before, 
since different frequencies are uncorrelated, separate matrix equations for each Fourier 

20 harmonic will be obtained. To accomplish this requires several steps. First, the AT 
independent elements of the matrix C at each frequency are listed as a column vector. 

Note that the order is not important. The outer product H(Q 0 )H\Q 0 ) is similarly 
reordered. The continuous integration is then approximated by a finite sum. There are, 
of course, an infinite number of ways to do this - the integral can be approximated using 
25 a trapezoidal or Simpson's rule, or expand D in a set of orthogonal functions. Perhaps 
the most natural expansion in this case is to write D as a Fourier series. 

+00 

D = Y a « eM A113 
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10 



15 



20 



Truncating this series at P terms, and substituting into equation A2.1 gives 

C = Ma Al. 14 

where a = (tf - y a + i) isaPxl column vector of directional Fourier amplitudes 

(the superscript "t" denotes transpose), Q is an N x 1 column vector, and 

M= r^GH(9)^ t (e)e" e 



A' (a) 



Al. 15 
Al. 16 



is a N x P matrix. From equation Al .6 we have 



11(9) =ae 



ik(r x cos9+r sin6) 



(F cos0 , F sin0 ,-/F') f 



where F(kz) = cosh (k(z + rf))/sinh(fo/), and F'(x) = dF/dx. We have used t = (r„rj = & 
-x'to denote the vector separation between range bins. 

Hence the expression for M in equation Al . 14 is seen to involve terms such as 

f 1 ^ 
cosG 

_ f +7t rfQ e inQ e ik ^ cos9 +r y sine > 

J-Jt 



1 

COS0 

sin0 
cos 20 
sin 20 



sin0 
cos 20 
sin 20 



Al. 17 



where r = ^rf + r^,(p =tan(r y /rj, and [...]„ denotes the function on the right-hand 

side of equation A1.19. 

The various terms above can be expressed in terms of the more primitive 

integral [e ipB ] n , which in turn can be evaluated in terms of Bessel functions as follows: 

[e ipe ]= V % dde^e^^^e^ Al. 18 

= 2n(iY +p e i(n ^J n , p (kr) Al. 19 

where J n (x) is a Bessel function of the first kind of degree n. 
Hence, 

[cos/ ? e] B =27i(0^^^[^J )I+/ ,(^)-(-)''^J (I+p (^)] A1.20 
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and 

[sin pQ] n = 2ji(/)" + ' ie* [e ip * J „ +p (kr) -(-)<> J n+p (kr) ] Al. 21 

These expressions permit M to be evaluated in terms of known functions at each 
frequency (and therefore each wavenumber). 

5 

A2. The Maximum-Likelihood Method (MLM) 

To form the MLM estimator, it is assumed that D can be expressed as a linear 
function of the data: 

Z)(0)-y(e) t Cy(e) A2. 1 

10 where the steering vector, y(9), is determined by minimizing the power in all directions, 
subject to the constraint that the signal associated with a plane wave is passed with unity 
gain. The solution for the weights can be obtained using a Lagrange multiplier and is 
given by 

v(B)= CH A2.2 

15 where C 1 is the inverse of C, which is computed using a singular- value decomposition. 
Note that y(0) is a function of both the array geometry/response and the data through its 
dependence on H and C. Then the MLM estimator can be written 

D(Q)= A2.3 
H^C' X H 

where the normalization constant, N D , is chosen so that D has unit area. The MLM 
20 estimator is "data-adaptive" is the sense that it uses the data to obtain filter weights y 
that for any given "look-direction" minimize the power coming from other directions. 
As a consequence, the expected value of the MLM estimate is the convolution of the 
true spectrum D with a window that depends on D. Since the window is a priori 
unknown, it is not possible to deconvolve it directly and another approach, such as the 
25 iterative maximum likelihood method (IMLM), must be employed. 

Although the off-diagonal elements of the data covariance matrix C are not 
biased, they are subject to sampling variability - their estimates converging to the true 
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cross-spectra as the degrees-of- freedom N c tends to infinity. This uncertainty translates 
into a mean-square error in the magnitude of D (Capon and Goodman, 1970) as 

var(Z)) = E[Df A2. 4 

N c -L + l 

where var(.) and £[..] denote the variance and expected value respectively, N c is the 
5 number of degrees-of-freedom in the estimate of C, and L is the number of range cells 
used (C has dimension L x Z). In addition to this sampling error, the diagonal elements 
of C are biased by the Doppler phase noise. The expected magnitude of this bias 
depends weakly on range through the acoustic signal-to-noise ratio, but to first order is 
roughly isotropic and independent of depth, and hence adds a constant to the diagonal 
10 elements of C. Intuitively, it is expected that this contribution will raise the overall 
level of the spectrum, and will tend to appear more strongly where the array sidelobes 
are largest. 

A3. The Iterative Maximum-Likelihood Method (IMLM) 

15 The IMLM estimator was proposed by Pawka (1983) for the analysis of 

pressure-gage-array observations of wave direction, and later applied to PUV 
observations by Oltman-Shay and Guza (1984). Krogstad et al. (1988) used the 
technique to obtain wave direction from velocity observations taken with a horizontal- 
beam sonar. 

20 The major difficulty with the MLM estimate is that it is not consistent, in the 

sense that the forward relation (2) applied to the MLM estimate D ML does not reproduce 
the observed covariance, C, and therefore the MLM method does not yield a true power 
spectral density. 

However, it is possible to define a heuristic iteration scheme that enforces 
25 consistency in an approximate way. Consider the nonlinear operator, N, defined by 

N:D forward > £ inverse ) jy^) A3> \ 

The problem is to find a spectrum D that satisfies 

N(D) = D ML , A3. 2 
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where D MLM denotes the original MLM estimate based on the observed covariance 
matrix, C. 

Krogstad et al. propose the following iteration scheme for the approximate 
solution of equation 1 1 

5 D n+1 =D„+X[D UL -N(D tt )] A3. 3 

where D 0 = D MLM , and X is a relaxation parameter that can be freely chosen. In practice, 
Krogstad et al. find reasonable convergence after 4-6 iterations for X's slightly greater 
than unity (typically in the range 1.1 - 1.4). Note that once the initial estimate has been 
computed, the iteration can proceed without further reference to the data covariance, C. 
10 In general, the IMLM estimate more closely follows the input spectrum than the 

MLM estimate, and the widths, especially in the region of the dominant peak, are also in 
better agreement. 

A4. The Eigenvector (EV) Technique 

15 The EV method (Johnson, 1982) is another high-resolution bearing estimator. 

The basic idea is to decompose the data matrix into a signal and an orthogonal residual, 
which we identify as the noise subspace. To formulate the problem more specifically, it 
assumed that m range cells are being used, and model the array output as 

d(t) = As(t) + n(t) A4. 1 

20 Here d is an mxl vector of observations, s is a pxl vector of (unknown) random signals, 
n is an mxl vector of Gaussian noise (assumed to be independent of the signals), and A 
is an mxp matrix characterizing the array response. Taking the Fourier transform of this 
equation and forming the cross-spectral matrix gives 

C = ASA t + N A4.2 
25 where S and TV are respectively the signal and noise covariances, and all quantities 
depend on frequency. 

The first considered is where the columns of A are linearly independent and 
there are fewer signals than data -e.g. p < m. The EV algorithm estimates the signal 



-33- 



directions by finding the orthogonal noise subspace. Operationally, the separation is 
accomplished by an eigenvector decomposition of the data matrix C. Since the 
eigenvectors are orthogonal and p < m, then m - p eigenvectors will be orthogonal to the 
signal subspace. We denote these eigenvectors and the corresponding eigenvalues by 

5 {ft j } and {Xj } respectively, where j = 1, ,«-/?. Then to the extent that the array 

response function H is a good representation of the signal, the projection of H on {hj } 
will vanish for angles corresponding to incident waves (e.g. signals) and will be non- 
zero otherwise. Similarly, any weighted sum of projections (or their magnitude- 
squared) will vanish, and the reciprocal will display sharp peaks at signal directions. 
10 Taking equal weights gives the MUSIC algorithm (Schmidt, 1986). The EV estimator 
uses the inverse of the noise eigenvalues as weights, and can be written as 

D = ^ = N ° A4.3 

Z-f y=i J I J 1 

where C N is the restriction of C to the noise subspace, and N D is a normalization 
constant. One might expect that this is a somewhat better estimator since it emphasizes 
15 eigenvectors that are more "noise-like" (e.g. those having smaller eigenvalues). 

As an aside, we note that in the limit that the number of signals tend to zero, the 
EV estimator in equation A4.3 becomes identical with the MLM estimator given in 
equation 19. 

In concluding this section, it is emphasized the fact that the EV algorithm is a 
20 bearing estimator. Its attraction lies in its simplicity and computational speed - to 
estimate the directional power density it must still be followed by the iterative IEV. 
Simulation results demonstrate that the EV can provide a useful estimate of the 
dominant wave directions when a lower bound on the noise dimension can be obtained. 
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WHAT IS CLAIMED IS: 

1 . A system for measuring the directional spectrum of waves in a fluid 
medium having a substantially planar surface, comprising: 

a sonar system having a plurality of transducers for generating respective acoustic 
5 beams and receiving echoes from one or more range cells located substantially within the 
beams, the beams being inclined at a non-zero angle with respect to the surface of the fluid 
medium; and 

a computer program executed by a processor for calculating the directional 
spectrum associated with the waves from the received echoes, wherein the computer 
1 0 program further utilizes a sensitivity vector as part of the calculation of the directional 
spectrum. 

2. The system of Claim 1 , wherein the sonar system comprises a broadband 
acoustic Doppler current profiler (ADCP). 

3. The system of Claim 1 , wherein the received echoes are related to the 
1 5 current velocity within the range cells. 

4. The system of Claim 1, wherein the transducers are arranged in the Janus 
configuration. 

5. The system of Claim 1, wherein the transducers are in a phased array 
configuration. 

20 6. The system of Claim 1 , wherein the calculation of the directional 

spectrum includes: 

calculating a non-directional wave height spectrum; 
calculating a cross-spectral matrix; 

calculating the directional spectra at each observed frequency; and 
25 calculating the dimensional directional spectrum form the non-directional 

wave height spectrum, the cross-spectral matrix, the directional spectra, and the 
sensitivity vector. 

7. The system of Claim 1 , wherein the fluid medium is at least in part water. 



-35- 



8. A Doppler sonar system for measuring the directional spectrum of waves 
in a fluid medium having a substantially planar surface, comprising: 

a signal generator for generating a transmitted signal; 

a plurality of transmitting transducers operatively connected to the signal 
generator, the transducers generating acoustic beams based on the transmitted 
signal and projecting the beams into the fluid medium, the acoustic beams further 
having an angular relationship to the surface of the fluid medium; 

a plurality of receiving transducers for receiving samples from one or more 
range cells located substantially within the acoustic beams, and producing a 
received signal relating to the samples; 

a signal processor capable of processing the received signal; and 

a computer program, executed at least in part on the signal processor, for 
calculating the directional spectrum associated with the waves based on the 
received signal, wherein the computer program further utilizes a sensitivity vector 
as part of the calculation of the directional spectrum. 

9. The sonar system of Claim 8, wherein the transmitting and receiving 
transducers are embodied within at least one unitary transducer element. 

10. The sonar system of Claim 8, wherein the transducers are arranged in the 
Janus configuration. 

11. A system for measuring the directional spectrum of waves in a fluid 
medium, comprising; 

a signal generator generating signals associated with one or more acoustic 

pulses; 

a plurality of transducers transmitting the acoustic pulses into the fluid 
medium in respective acoustic beams, the transducers further receiving echoes 
generated by the acoustic pulses from one or more range cells located within the 
acoustic beams; and 

a signal processor receiving signals indicative of the received echoes, and 
calculating a sensitivity vector associated with the first and second transducers. 
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1 2. The system of Claim 1 1 , wherein the acoustic beams are inclined at a non- 
zero angle with respect to the surface of the fluid medium. 

13. The system of Claim 1 2, wherein the acoustic beams project downward 
from the surface of the fluid medium. 

5 14. A method of calculating the directional spectrum of a wave in a fluid 

medium utilizing a plurality of acoustic transducers, comprising : 

generating a plurality of acoustic beams from the transducers, the beams 
having range cells located at least partly therein; 

receiving echoes produced within an array consisting of two or more of 
10 the range cells; 

processing signals indicative of the received echoes; 
generating a sensitivity vector, the vector being based on the signals and 
the geometry of the array of range cells; and 

estimating the directional spectrum of the wave based on the signals and 
1 5 the sensitivity vector. 

15. The method of Claim 14, wherein the processing of the signals indicative 
of the received echoes includes: 

decoding the raw data associated with the received echoes; 
calculating the surface height for each acoustic beam; 
20 calculating a mean value of current in the fluid medium; and 

numerically inverting a linear dispersion relation, wherein the linear 
dispersion relation relates wave frequency, water depth, and wavenumber. 

16. The method of Claim 14, wherein the generating of the sensitivity vector 
includes: 

25 selecting at least one range cell from at least two of the acoustic beams; 

calculating a plurality of velocity components for each of the selected 
range cells; 

calculating a plurality of Fourier coefficients associated with each of the 
velocity components; and 



-37- 



calculating a sensitivity vector from the plurality of Fourier coefficients. 

1 7. The method of Claim 1 4, wherein the estimating of the directional 
spectrum includes: 

calculating the wave height spectrum S H ; 
5 calculating the cross-spectral matrix C; 

calculating a directional spectrum at each observed frequency; and 
constructing the estimate of the dimensional wave directional spectrum 
from the directional spectra. 

1 8. The method of Claim 14, wherein the sensitivity vector includes elements 
1 0 corresponding to surface height and pressure within the fluid medium. 

1 9. The method of Claim 1 4, wherein the estimating the directional spectrum 
comprises maximum likelihood processing of the signals indicative of the received 
echoes. 

20. The method of Claim 1 9, wherein the estimates of the directional spectrum 
1 5 further comprises iterative maximum likelihood method (EMLM) processing. 

21. A method of measuring the directional spectrum of waves in a fluid 
medium using a sonar system having an upward or downward looking transducer 
configuration, comprising: 

generating one or more acoustic beams from the transducer configuration; 
20 measuring the current velocities within one or more range cells of the 

acoustic beams; 

forming a sensitivity vector related to the transducer configuration using a 
linear wave model; and 

forming a wave directional spectrum matrix using the measured current 
25 velocities and the sensitivity vector. 

22. A computer program used with a sonar system for calculating the two- 
dimensional directional spectrum of a wave in a fluid medium from at least one set of 
received echoes, comprising: 
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a data processor for processing the signals representative of the received 

echoes; 

means for calculating a non-directional wave height spectrum from the 
signals processed by the data processor; 
5 means for generating a cross-spectral matrix from the signals; 

a sensitivity vector generator, wherein a sensitivity vector is generated 
relating to the configuration of the transducers of the sonar system; and 
a two-dimensional directional spectrum estimator, the estimator 
calculating the directional spectrum from the non-directional wave height 
1 0 spectrum, the cross-spectral matrix, and the sensitivity vector. 

23 . A method of measuring the directional spectrum of waves in a fluid 
medium using a sonar system having an upward or downward looking transducer 
configuration, comprising: 

generating a plurality of acoustic beams from the transducer configuration; 
1 5 measuring current velocities within one or more range cells of the acoustic 

beams; 

calculating values of wave frequency and wave number magnitude 
according to a linear wave dispersion relation; 

forming a wave height spectrum matrix using the measured current 

20 velocities 

forming a cross-spectral matrix as a function of wave frequency; 
generating a sensitivity vector related to the transducer configuration of the 
sonar system; 

using a maximum likelihood method to generate directional spectra for 
25 each observed frequency from the cross-spectral matrix and sensitivity vector; and 

constructing a complete two-dimensional wave directional spectrum from 
the directional spectra for each observed frequency. 

24. The method of Claim 23, wherein the construction of the complete two- 
dimensional wave directional spectrum further includes: 
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renormalizing the observed frequency-dependent wave directional 
spectrum; 

determining the presence of a non-zero current; and 
in the event of non-zero current, remapping the frequency and applying a 
5 Jacobian operator. 

25. The method of Claim 23 , wherein the maximum likelihood method used in 
generating the directional spectra further comprises an iterative maximum likelihood 
method. 

26. The method of Claim 23, wherein the sensitivity vector includes elements 
1 0 corresponding to surface height and pressure within the fluid medium. 

27. The system of Claim 1 , wherein the processor is independent from the 
sonar system. 

28. The system of Claim 1 , wherein the processor is included in the sonar 

system. 

15 29. The system of Claim 1 , wherein the processor comprises a signal 

processor. 



20 
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SYSTEM AND METHOD FOR MEASURING 



WAVE DIRECTIONAL SPECTRUM AND WAVE HEIGHT 



Abstract <rf the Disclosure 
5 A system and method for measuring the directional spectrum of one or more 

waves in a fluid medium using a multi-beam sonar system. In an exemplary embodiment, 
range cells located within a plurality of acoustic beams are sampled to provide current 
velocity data. Optionally, wave surface height and pressure data is obtained as well. This 
velocity, wave height, and pressure data is Fourier-transformed by one or more signal 

10 processors within the system, and a surface height spectrum produced. A cross-spectral 
coefficient matrix at each observed frequency is also generated from this data. A 
sensitivity vector specifically related to the ADCP ! s transducer array geometry is used in 
conjunction with maximum likelihood method (MLM), iterative maximum likelihood 
method (IMLM), or other similar methods to solve a the wave equation at each frequency 

15 and produce a frequency-specific wave directional spectrum. Ultimately, the frequency- 
specific spectra are combined to construct a complete two-dimensional wave directional 
spectrum. The system is also capable of measuring current profile as a function of depth 
in conjunction with wave direction and wave height. 

20 S:\DOCS\RFG\RFG-1172.DOC 
080498 
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